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Abstract 



A strongly focused laser beam can be used to trap, manipulate and exert torque on a microparticle. The 
torque is the result of transfer of angular momentum by scattering of the laser beam. The laser could be used to 
drive a rotor, impeller, cog wheel or some other microdevice of a few microns in size, perhaps fabricated from a 
t-^ ■ birefringent material. We review our methods of computationally simulating the torque and force imparted by 

a laser beam. We introduce a method of hybridizing the T-matrix with the Finite Difference Frequency Domain 
(FDFD) method to allow the modelling of materials that are anisotropic and inhomogeneous, and structures 



that have complex shapes. The high degree of symmetry of a microrotor, such as discrete or continuous 
rotational symmetry, can be exploited to reduce computational time and memory requirements by orders of 
magnitude. This is achieved by performing calculations for only a given segment or plane that is repeated 
across the whole structure. This can be demonstrated by modelling the optical trapping and rotation of a 
'"5 ! cube. 

ft ; 

1 Introduction 

> ■ 

The T-matrix method |lj is commonly used to calculate properties of light scattering from axisymmetric 
^vq ' mesoscale (jA — 5A) particles that are homogeneous and isotropic [2J |3|. Using the T-matrix, the optical force 
[~^. . and torque imparted on the particle by the incident beam can be calculated |3|. The T-matrix is independent of 
the incident field and only dependent of the properties (size, shape, orientation, permittivity) of the particle. If 
the incident fields change, the T-matrix need not be recalculated. 

We extend this method to model particles that are inhomogenous, anisotropic and have complex geomet- 



rical shapes by combining the T-matrix method with the Finite Difference Frequency Domain (FDFD) method. 
• i-H In the FDFD method, we discretize the computational region into a grid with sufficiently small grid size. The 
inclusion of FDFD equations in the algorithm is computationally intensive. To optimize computational time 
and memory usage, we consider the rotational symmetry of the system. If the particle is rotationally symmet- 
ric about an axis, the system could be reduced to a 2D problem (figures la and lb). By choosing a cylindrical 
coordinate system, the section of interest can be treated in 2D rectangular (r,z) coordinates which leads to 
compatibility with the FDFD cell (figure lc). If the particle has nth-order discrete rotational symmetry, typical 
of a microrotor, savings in computational time and memory could still be achieved by performing calculations 
for only one repeated segment. For example, we modelling the optical trapping of a cube, exploiting the 4th 
order rotational symmetry and xy-plane mirror symmetry, to reduce the time required to calculate the T-matrix 
from 30 hours to 20 minutes. 



2 FDFD equations 

The rotationally symmetric FDFD equations were derived by expanding Maxwell curl and divergence equa- 
tions 1 4 1 in cylindrical coordinates 



VxA= ^ f+ ^ H -^"^ t (1) 
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Figure 1: a) Cylindrical coordinate system. Spheroidal particle enclosed in a cylindrical volume, b) Rectangular 
computational grid with regions (1) and (2), inside and outside the particle respectively, c) FDFD (Yee) cell. 
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where A represents the electric field, E, or the magnetic field, H. The evolution of both fields can be expressed 
as dA/dt = — icuA. It is sufficient, as will be seen below, to consider a field with azimuthal variation exp(im$); 
the variation of the field with respect to <p would be dA/d(p = imA. Substituting the time evolution and 
4> harmonic equations into Maxwell's equations we obtain 6 curl equations and 2 divergence equations for 
electric and magnetic fields. As an example, the curl equation for E r is 
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Using the Yee cell 1 5 1 in figure lc for discretization, the FDFD curl equation for E r can be written as 
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The other equations can be discretised similarly. 



3 Hybridizing the T-matrix method with FDFD 

The T-matrix is an operator (T) which acts on the coefficients of the incoming field to produce the coefficients 
of the outgoing field 

p = Tfl (5) 

where a represents the vector made up of the coefficients (a nm and b nm ) of the incoming field and p represents 
the vector of the coefficients (p nm and q nm ) of the outgoing field. The electric fields (and similarly for magnetic 
fields) can be expanded in terms of incoming and outgoing Vector Spherical Wave Functions (VSWFs) 

oo n 

^in = 7. 7. a nm MjJ (k out r) + btm^nm (hutf), (6) 
n=l m=—n 

Eout = L PnmM^ikoutr) +CjnmNliJ(koutr). (7) 

n=\ m=—n 

where k ut is the wave vector outside the particle, and M and N are vector spherical wave functions (VSWFs) 
defined in Q. Naturally, we cannot take the sums to infinity but rather taken to N max which is based on 
criteria defined in 1 2 1 . In our model, we would have a dielectric region within the computational grid that 
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would interact with the incoming and outgoing fields. So, in coupling the electric field E(r) from the FDFD 
solutions with the VSWFs for the TE incident modes we obtain 

^2U r ) + £ PnmM { „%) + q nm N^(r) = E(r), (8) 
n=l 

where n' is the incident mode. Similarly for the TM modes, 

tiSl ( r ) + E P«mN^ (r) + q„ m M^2 (r) = E(r). (9) 

n=l 

Due to the rotational symmetry, there is no coupling to other azimuthal modes (i.e. only one value of m' 
appears). Therefore, all fields share an azimuthal dependence of exp(im$). Equation (8) or (9) connects the 
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Figure 2: a) Coefficient Matrix, Unknowns and Constants, b) The T-matrix is made up of four blocks that result 
from the coupling between the TE and TM modes. The inset is a section of the T-matrix in detail where the 
allowable coupling between m and n modes are shown as black pixels. 



VSWF description of the external fields and the FDFD grid. The VSWF and FDFD equations form an overde- 
termined linear system (figure 2a) and can be solved using a standard numerical library. The FDFD equations 
are inserted in the Coefficient matrix (figure 2a) first followed by VSWF equations. Last, the z-axis boundary 
equations are inserted in the Coefficient matrix. Generally, the field is zero at the z-axis except for the modes 
m = ±1, in which case the first derivatives of the fields are zero. Cycling through all incident modes, the 
solutions for p nm and q nm are solved given one incident mode at a time and their values are inserted into the 
T-matrix column representing coupling between the m and n modes (figure 2b). 



4 Discussion 

The micromachines of interest may or may not have xy-plane mirror symmetry but they will typically have 
nth-order rotational symmetry. Nonetheless, as with the cube we had modelled, the rotational symmetry 
can be exploited to reduce the calculation time by orders of magnitude. Conventional T-matrix methods are 
limited in their application to modelling homogeneous and isotropic materials, with shapes that are close 
to spheroidal. The FDFD hybridization extends the modelling capability to include nth-order rotationally 
symmetric micro-machines with complex shapes made from materials that are inhomogenous and anisotropic 
e.g. birefringent crystals. 

The Matlab script for solving the matrices in figure 2a was tested on a PC with a 32-bit single 3GHz CPU 
and 1Gb of RAM. We performed the calculation simulating a 3000 nm radius cylinder with grid sizes from 
1000nm-250nm. Extrapolating from the natural log scale plot (figure 3), we estimated that it would take 13.6 
hours and 165.9 hours (7 days) to perform the calculation given 100 nm and 50 nm grid spacing respectively. 

While the foregoing is directed at modelling rotational symmetric particles, we intend to model more com- 
plex particles buy using a 3D FDFD grid or the Discrete Dipole Approximation (DDA) method 1 1 [ coupling 
with the VSWFs. 
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Figure 3: Log gridsize in wavelength units versus log time (sees). 
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